Extended Kramers-Moyal analysis applied to optical trapping 
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The Kramers-Moyal analysis is a well established approach to analyze stochastic time series from 
complex systems. If the sampling interval of a measured time series is too low, systematic errors 
occur in the analysis results. These errors are labeled as finite time effects in the literature. In the 
present article, we present some new insights about these effects and discuss the limitations of a 
previously published method to estimate Kramers-Moyal coefficients at the presence of finite time 
effects. To increase the reliability of this method and to avoid misinterpretations, we extend it by 
the computation of error estimates for estimated parameters using a Monte Carlo error propagation 
technique. Finally, the extended method is applied to a data set of an optical trapping experiment 
yielding estimations of the forces acting on a Brownian particle trapped by optical tweezers. We 
find an increased Markov-Einstein time scale of the order of the relaxation time of the process which 
can be traced back to memory effects caused by the interaction of the particle and the fluid. Above 
the Markov-Einstein time scale, the process can be very well described by the classical overdamped 
Markov model for Brownian motion. 

PACS numbers: 05.10.Gg, 05.45.Tp, 05.40.Jc, 87.80.Cc 



I. INTRODUCTION 



Many real world stochastic processes q(t) can be mod- 
eled by stochastic differential equations of Langevin type 



respectively. Given the transition probability densities 
p q (x' ,t + r\x,t) of a general Markov process q(t), the nth 
KM coefficient can be defined as 



D {n) (x) = lim ilWfa 

r^O nW 



(4) 



q = DW(q) + yj2DW(q) 



r. 



(1) 



where D^ x \q) is called the drift coefficient that describes 
deterministic influences on the dynamics and D^ 2 \q) is 
the diffusion coefficient that describes the state depen- 
dent amplitude of a fast fluctuating force r(t) with zero 
mean, adopting Ito's interpretation. A frequently applied 
idealization is to assume that r(t) is S correlated in time 
and Gaussian distributed. In this case the process fulfills 
the Markov property and can be completely characterized 
by the corresponding Fokker-Planck equation (FPE) 



where Mr (x) is called the nth conditional moment 
given by 



M^\x) = {{q{t + T)-q(t)T)\ q{t)=x 



(x' - x) n p q (x', t + t\x, t)dx' 



(5) 



d_ 

at 



fg(x,t) = L(x)f g (x,t) 



with the Fokker-Planck operator 
L(x) = 



(2) 



(3) 



that describes the temporal evolution of the probability 
density function (PDF) f q (x,t) = (S(x - q(t))) of the 
process q(t). The Fokker-Planck operator L(x) depends 
on the drift and diffusion coefficients that determine the 
corresponding Langevin equation (Q~]). 

Drift and diffusion coefficients are also referred to as 
the first and second Kramers-Moyal (KM) coefficients, 
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Since the transition PDFs can be estimated from mea- 
sured time series data of a process of interest, it is in 
principle possible to set up a model in terms of a Fokker- 
Planck or Langevin equation by data analysis. Since 
Friedrich and Peinke applied this nowadays called KM 
analysis to the investigation of the turbulent cascade 
[J, Q in 1997, it has developed to a rapidly growing field 
of research with many applications in the natural sciences 
and beyond (see Q and references therein). 

One major problem connected to the KM analysis is 
the limit of the time increment r to zero that has to be 
performed in the determination of the KM coefficients. 
Without an appropriate limiting procedure, estimated 
KM coefficients can deviate significantly from the true 
coefficients if the minimal available r, i. e. the sampling 
interval of the measured time series, is too large. These 
errors are referred to as finite time effects and were the 
subject of many publications during the last years (e. g. 
[7t-lll1] ) . Without an adequate understanding of finite 
time effects, KM analysis involves the risk of dramatic 
misinterpretations of the achieved results. 

In a recent publication we have developed a 

method that allows for a correct KM analysis when the 
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sampling interval is large. In the present article we 
demonstrate that also this method fails if the sampling 
interval is so large that the information about the true 
KM coefficients is (almost) lost. In order to lower the 
risk of misinterpretations in those cases, we extend the 
method by the computation of error estimates for the ob- 
tained model parameters. The problem of error estimates 
for model parameters in the KM analysis context is also 
addressed in a very recent publication by Kleinhans [TJ[ 
in a Bayesian framework. In contrast to our approach, 
this method is only designed for data sets with a suffi- 
ciently large sampling rate. 

We present an application of our method to an opti- 
cal experiment that yields real-world stochastic data via 
trajectories of a Brownian particle trapped by optical 
tweezers. We apply our method to these trajectories in 
order to yield estimations for the forces induced by the 
optical tweezers and to analyze the spatial distribution 
of the temperature. 

The outline of the article is as follows. In Sec. |TT] we 
review the problem of finite time effects and discuss an 
example in order to provide an intuitive understanding. 
Afterwards we review the method of Ref. [l2j and present 
important extensions. Finally, Sec. IIVI describes the ap- 
plication to our experiment. The last section is dedicated 
to some concluding remarks. 



II. FINITE TIME EFFECTS 

At first we introduce the finite time KM coefficients 

flW(i) = iMW(i), (6) 
n'.T 

where (x) are the conditional moments defined in 
Eq. ©. Apparently, lim r ^ D { T n) (x) = D^ n \x). In 
order to understand finite time effects and to possibly 
correct for them, it is important to know, how the finite 
time coefficients change with the time increment r given 
the true coefficients. 

There are several possibilities to compute finite time 
coefficients. One possibility is by solving the FPE. The 
transition PDF p q (x, to + t\xq, t ) that occurs in Eq. §5§ 
is the solution to the FPE @ at time t = to + t with the 
initial condition 



(7) 



Unfortunately this is in most cases impossible to do an- 
alytically. Also for a numerical solution of the FPE, the 
initial condition in form of a Dirac S distribution will 
cause problems. 

Another possibility is to numerically integrate the cor- 
responding Langevin equation and to estimate the finite 
time coefficients from the generated time series. Based 
on this approach, Kleinhans et al. developed a maxi- 
mum likelihood approach to estimate KM coefficients for 
processes with finite sampling rates [IJ], [U| • 



A third possibility is the adjoint operator approach. 
In Ref. [8j, Friedrich et al. mention that the conditional 
moments can be expressed as a Taylor series in r which 
reads 



M^\x) = 



E 



fc! 



(x' - xf 



(8) 



This series expansion contains the adjoint Fokkcr-Planck 
operator 

LHx')^D^(x')^+D^(x')^. (9) 

With the adjoint operator method it is not only possible 
to obtain a series expansion. As Lade showed in Ref. 
[9J, the conditional moments can also be computed by 
solving the partial differential equation 



dW n , x (x',t) 
dt 



= D(x')W n<x {x l ,t). 



(10) 



(11) 



(12) 



Eq. ([10]), also called the adjoint FPE (AFPE), can easily 
be solved analytically as long as the drift is linear. Oth- 
erwise numerical solutions can be obtained with standard 
finite difference schemes. 

For an Ornstein-Uhlenbeck (OU) process, 



with initial condition 

W n , x (x',0) = (x'-x) n . 
Then the conditional moments are given by 
M<p(x) =W n>x (x' = x,t = r) . 



D^\x) = -jx. 
D ( - 2 \x)=a, 



(13a) 
(13b) 



the analytic expressions for the finite time coefficients 
read 



1 



W0 = £ 



-e" 7r ) 
fl-e^ 



T + -(i- 

7 



-2jt\ 



(14a) 



(14b) 



We have plotted the finite time coefficients together 
with the conditional moments in Fig. [Tj (drift) and Fig. 
[2] (diffusion) for 7 = a = 1. By taking a logarithmic scal- 
ing for the r axis, one can clearly identify two limiting 
cases separated by the relaxation time, which is tr = 1 
in this case. For r <C tr, the finite time KM coefficients 
have converged to the true KM coefficients, whereas the 
conditional moments tend to zero. For r 3> tr, the con- 
ditional moments become stationary with respect to r, 
whereas the corresponding KM coefficients vanish. 

The latter case is the limit of statistical independence 
that has been discussed by Anteneodo and Queiros in 
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FIG. 1. Conditional moment (top) and finite time drift (bot- 
tom) for an OU process with D* 1 ' = —x and — 1. The 
r axes are in logarithmic scaling. 

Ref. [lOj. For r tr, the realization a; of a process 
q at time t becomes statistically independent from the 
realization x' at time t + r. Therefore, 

Pq (x',t + T\x,t)ttf q (x',t + T). (15) 

For a stationary process, f q (x', t+r) is just the stationary 
PDF f q {x). As a consequence, the conditional moments 
yield 

MW TR (x)~{q)-x, (16a) 
M™ TR (x)~ (q 2 ) -2{q)x + x 2 . (16b) 
and the finite time KM coefficients read 

D% TR {x)^-[{q)-x] , (17a) 

T 

D?l TR (x)^^-[(q 2 )-2(q)x + X 2 ] . (17b) 

That means that every Langevin process, independent 
of the true KM coefficients, appears as a process with lin- 
ear drift and quadratic diffusion if the sampling interval 
of the data set is large compared to the relaxation time of 
the process. Therefore, as a preinvestigation before a KM 
analysis, it is necessary to investigate the autocorrelation 
function of a process in order to estimate the relaxation 
time, and to estimate the conditional moments for var- 
ious available r and present it as a plot like Figures Q] 




FIG. 2. Conditional moment (top) and finite time diffusion 
(bottom) for the same process as in Fig. Q] 



and [21 Furthermore, Ref. pj| provides some checks in 
order to validate models with linear drift and quadratic 
diffusion. 

We can conclude that in the regime r <C tr, finite 
time effects are small and low order corrections should 
suffice or can even be neglected. In the regime r 3> tr, 
all information about the true dynamics is lost and a 
KM analysis is no longer possible. In between those two 
regimes, if t rj tr, a KM analysis is possible with use 
of our method presented in Ref. (H| , as is demonstrated 
in several examples. Of course, at some point depending 
on the available amount of data, the results obtained by 
this method will also fail when r approaches the limiting 
case of statistical independence. The point at which it 
fails will depend on the size of the available data set. 
The lower statistical uncertainties are in the estimation 
of the conditional moments, the more accurate one can 
extrapolate to r — 0. To quantify this possible accuracy, 
we extend our method by the computation of errors in 
the estimated parameters with which drift and diffusion 
coefficients are parametrized. This is done by a Monte 
Carlo error propagation technique that is often used in 
nonlinear optimization problems [16j . 

In the next section, we first summarize the general 
method and present some minor improvements. After- 
wards the Monte Carlo error propagation technique will 
be introduced in some more detail. We also address the 
limitations of the method due to finite amount of data 
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and large sampling intervals. 



the conditional moments. In 11211 we used 



III. KM ANALYSIS AT LOW SAMPLING 
RATES 

A. Review of the method and minor improvements 

The optimization procedure introduced in [l2[ is based 
on an inversion of Lade's method [9J to compute finite 
time coefficients with help of the AFPE, which was al- 
ready mentioned in the foregoing section. The basic steps 
that have to be performed are the following. 

As a first step, one has to obtain estimations of the 
conditional moments M T ' (x) for several values of x and 
t from the time series data. For this we recommend a 
kernel based regression with use of the Nadaraya- Watson 
estimator. The bandwidth was in [ljj selected according 
to Silverman's rule of thumb [I?} as 



h = \m&N- 1/5 . 



(18) 



where a is the standard deviation of the time series data. 
We found that this bandwidth is in many cases a bit too 
large leading to underestimated slopes of straight lines or 
curvatures of parabolas. More sophisticated data driven 
bandwidth selectors such as cross validation techniques 
or the method described in [1 81 ] are computationally ex- 
pensive and did not yield robust improvements in our 
tested cases. We obtained the best results by simply re- 
ducing the factor 1.06 in ([TS]) to 0.8. 

The next step is to find a suitable parametrization for 
the drift and diffusion coefficients. In most cases a poly- 
nomial ansatz will be appropriate. Also a representation 
in form of spline curves has been successfully tested in 
[l2| , where the interpolation points serve as optimization 
parameters. In either case, the coefficients are repre- 
sented as D( 12 \x, a) with a set of optimization param- 
eters a. If an analytic expression has been chosen for 
the parametrization, one can try to obtain an analytic 
expression for the corresponding conditional moments 

(1 2) 

Mr ' (x,a) via Lade's method. Otherwise, they have 
to be calculated by numerical integration of the AFPE. 
The optimal set of parameters must then be obtained by 
minimizing a distance measure between the estimated 
conditional moments M T ' (x) and the computed ones. 
In [l2} we have chosen the weighted least squares distance 



M N 

1=1 j=i 



(19) 



,(2) 





M Ti '(xj)j 








M T f{x 3 )-[ 




ELi \k 







(20a) 



(20b) 



as an error estimate, where T is the number of data points 
of the time series data, K(») is the kernel function, and h 
is the selected bandwidth. The error estimate ([20)) does 
not take into account the influence of the bandwidth. 
According to [13], Eqs. (|20|) have to be corrected by a 
factor of 



\K\ 



where H-SfH^ is the L2 norm of the kernel, thus 



\K\\ 2 = 



K 2 (x)dx 



For the Epanechnikov kernel 



K(x) 




5-x 2 



if 



> 5 



(21) 



(22) 



(23) 



(1 2) 

where a) „■ ' are the statistical errors in the estimation of 



which we use, ||if||2 ~ 0.518. However, regarding the 
optimization, the correction factor (|21[) has no influence, 
since the error estimate has only the purpose of a rela- 
tive weighting of the estimated conditional moments for 
different values of x and r. It only becomes relevant in 
connection with the Monte Carlo error propagation de- 
scribed in Sec. IIII B[ 

The minimum of (|19[) is in [l2[ determined by a trust 
region algorithm, but other optimization methods should 
work as well. 



B. Monte Carlo error propagation 

The basic idea of the Monte Carlo error propagation 
(MCEP) approach is as follows: After estimating the 
model parameters of a specific model from a noisy data 
set (backward problem), the model is used to produce 
an unnoisy data set (forward problem) . Then one gener- 
ates an ensemble of pseudo data sets by perturbing the 
unnoisy data set with different realizations of noise. For 
each data set again, the model parameters are estimated. 
Eventually, one can compute the standard deviation for 
each model parameter from this ensemble, which can be 
used as an uncertainty measure. 

In the present case, we do not regard the time series 
data as our noisy data set, but the estimated conditional 
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moments M T ' (x). In order to obtain reliable uncer- 
tainty estimates with the MCEP approach, it is crucial 
that the artificial perturbations of the unnoisy data set 
are of the same amplitude as the errors of the original 
data set. Therefore one needs a reliable absolute error 
estimate for the conditional moments and the correction 
factor (f2Tj) becomes important. 

To generate the ensemble of pseudo data sets, one has 
to compute the conditional moments that correspond to 

the estimated model parameters by solving the corre- 
al 2) 

sponding AFPE. To each M Ti ' (Xj) we add an indepen- 
dent Gaussian distributed random number with standard 
deviation of expression ([20]) multiplied by the factor (j2"Tj) 
for each member of the data set ensemble. Then the 
model parameters are estimated again for each data set 
and the corresponding standard deviation is calculated 
for each model parameter. 

There is one technical problem connected to the MCEP 
approach. It can happen that the optimization routine 
does not find the absolute minimum of the least square 
potential for some pseudo data sets of the ensemble. Al- 
though the trust region algorithm is a very robust one, 
our observation is that it sometimes gets stuck in local 
minima close to the initial condition. If the initial con- 
dition is selected the same for all pseudo data sets, it 
can therefore happen that the computed standard devi- 
ation is significantly underestimated. To overcome this 
problem, we choose random initial conditions that lie in a 
region around the estimated parameters from the original 
data set. This in return bears the danger that the opti- 
mization routine will not find the way back to the true 
minimum if the initial condition lies too far away. This 
can lead to overestimated standard deviations. There- 
fore, only optimization results are accepted where the 
final residual is below a specific threshold. Otherwise 
the optimization is repeated with another initial condi- 
tion. The threshold is selected three times the residual 
of the result from the original data set. Nevertheless 
it can happen that some outliers lead to a significantly 
overestimated standard deviation, especially in cases of 
very large sampling intervals. The easiest way to over- 
come this problem is to detect these outliers by eye and 
remove them from the ensemble. 



TABLE I. Comparison between the standard deviations cje of 
estimated parameters in the ensemble of ten OU processes and 
the mean errors omcep computed by the MCEP approach. 



Parameter 


average 


<JE 


CMCEP 


7 


1.0014 


0.0035 


0.0033 


Q 


1.009 


0.012 


0.014 


ft 


-0.0006 


0.0063 


0.0116 



To test the MCEP approach, we generate an ensemble 
often synthetic time series of an OU process with — 
— x and — 1. The corresponding Langevin equation 
is integrated with the Euler-Maruyama scheme with a 
time increment At = 10~ 2 . Only every 100th data point 
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FIG. 3. Estimated parameters for an ensemble of ten OU 
processes according to Eqs. (|24|l with different noise realiza- 
tions. The error bars are computed via the MCEP method. 
The horizontal dashed lines indicate the true values of the 
parameters. The sampling interval r of the series is equal to 
the relaxation time tr = 1. 

is stored so that the sampling interval of time series is 
t = 1. Each time series consists of 10 6 data points. For 
the optimization we choose the parametrization 

D {1 \x) = -jx, (24a) 

D (2) {x) = a + ftx 2 , (24b) 

with the optimization parameters 7, a and ft. For each 
of the ten data sets, we estimate these parameters as 
well as the corresponding standard deviations with the 
MCEP approach. The results are depicted in Fig. [3] 
By eye, the estimated error bars seem reasonable. To 
make a quantitative statement, we compute the standard 
deviations of the estimated parameters from the ensemble 
and compare them to the mean standard deviations of the 
MCEP approach. The numbers are listed in Table fl] The 
errors of the parameters 7 and a fit very well while the 
error of ft is a bit overestimated. 

However, one should note that the errors obtained 
by the MCEP approach only cover uncertainties caused 
by statistical fluctuations due to the finite amount of 
data. The influence of other error sources such as mea- 
surement noise, deviations from the Markov property, 
non-stationarity of the time series or an inappropriate 
parametrization of drift and diffusion coefficients, just to 
name a few, are not considered. Therefore, the MCEP 
errors should be regarded as a lower bound for the true 
errors. Nevertheless they can decrease the danger of an 
overestimation of the significance of KM analysis results. 
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C. Limitations of the approach 

Now we use the MCEP approach to demonstrate the 
limitations of the KM analysis caused by finite time ef- 
fects and limited amount of data. To this end, we com- 
pute the error estimates for the parameters 7, a and j3 
for synthetic time series with different sampling intervals. 
The synthetic time series are generated in the same man- 
ner as described in section UlIBI All data sets consist of 
10 6 data points. 
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FIG. 5. Experimental setup: (D)M: (dichroic) mirror; C: 
CMOS camera; I: illumination; MI: microscope; MO: micro- 
scope objective; Tl and T2: telescope systems to achieve an 
optimal beam diameter; SLM: phase-only spatial light modu- 
lator, for the described experiment not connected to a power 
supply and thus acting as a mirror; SP: sample plane, where 
20 uL of the used suspension was fixed between a microscope 
slide and a coverglass. 



FIG. 4. (Color online) Estimated MCEP errors for synthetic 
time series according to the model (|24l) with parameters 7 = 
a = 1 and j3 — 0, as a function of the sampling interval r. 

The result is shown in Fig. 01 One can see that the 
estimation of the diffusion parameters a and /3 breaks 
down around 2.5 relaxation times (tr — 1), while the 
drift parameter 7 can still be estimated up to about 5 
relaxation times. Above these values it is difficult to ob- 
tain robust error estimates with the MCEP approach be- 
cause the afore-mentioned outliers become very frequent. 
However, using the MCEP approach for sampling inter- 
vals above these values, one will clearly notice that the 
KM analysis is no longer feasible. 

Looking at Fig. [4j one further notices that the errors of 
a and j3 are strongly correlated, which is not surprising. 
That the absolute sizes of the errors are almost equal 
is only the case for the specific selection of parameters 
7 = a = 1 and f3 = 0. 



IV. APPLICATION TO AN OPTICAL 
TRAPPING EXPERIMENT 

A. Experimental setup 

Trapping of small dielectric objects in the nano- and 
microscale can be realized by focusing a single laser beam 
so that it creates a gradient force near the focal region 
that is able to hold these particles, and even levitate them 
against gravity jioj . Our trapping experiment is based 
on such a single beam gradient laser tweezers system (see 
FIG. [5]) as in [2(| ■ A near infrared laser (A = 1064 nm, 
maximum output power of 2.5 W, Smart Laser Systems) 



with a Gaussian beam distribution is coupled into an in- 
verted microscope (Ti Eclipse, Nikon) for trapping with 
an oil-immersion objective (CFI Apo TIRF lOOx) with a 
numerical aperture of NA = 1.49. This objective is also 
used for observation purposes and images the observation 
plane onto a high-speed CMOS sensor (MV2-D1280-640- 
CL-8, Photonfocus) with a frame rate of 488 Hz at full 
resolution. By reducing the size of the captured images 
to a region of interest of 80x64 pixels, the framerate 
could be increased to 3966 Hz, which corresponds to a 
sampling interval t s s» 0.25 ms. For the measurements 
uncoated polystyrene beads (Kisker) with a diameter of 
(1.002 ± 0.043) |xm are used. These microparticles are 
suspended in destilled water and trapped 5-10 urn above 
the cover glass. Due to technical reasons up to 150 videos 
consisting of 10 4 frames are captured over 15 minutes 
and analyzed afterwards. The particle's position is de- 
termined by a Matlab algorithm implementing a center 
of mass detection [U • 

The autocorrelation function of a trapped particle's 
x(y)-position decays with a relaxation time T x f y \ — 
GTrr/r / k x (y) where 77 is the viscosity of the surrounding 
fluid, r the radius of the particle and k x M the corre- 
sponding trap stiffness [23]. For a particle with a ra- 
dius of r = 0.5 |_im typical relaxation times are around 
10 ms \1mpN~ 1 /k x . Therefore the laser power is reduced 
to about 10 mW in the trapping plane which results in 
a trapping stiffness of k x — 21pN/|xm and thus a relax- 
ation time of about t x = 0.45 ms (compare Fig. [5]) in 
order to have sufficient temporal resolution to observe 
the particle's dynamics. 
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B. Modeling 

The equation of motion for a Brownian particle in an 
optical trap reads for one spatial dimension [25| 



m p x(t) = + fop (*(*)) + Fth(t) 



(25) 



where m p is the mass of the particle, F{ r is the friction 
force, F op is the force induced by the optical trap and 
.Fth denotes the thermal fluctuations. The generell form 
of the friction force reads 



j(t - t')x(t')dt' . 



(26) 



According to the fluctuation dissipation theorem [24( , the 
kernel 7(t) is connected to the correlation function of the 
thermal fluctuations 



{F ih (t)F th (t'))=k B T 7 {\t-t'\). 



(27) 



Assuming a laminar velocity profile around a spherical 
particle and no-slip boundary conditions, the kernel is 
given by 



j(t - t') = 2S(t - t')6TTrjr 



(28) 



where r\ is the dynamic viscosity of the fluid and r the 
radius of the particle. This yields the well-known Stokes' 
law 



F fr (t) = -Xx(t) 



(29) 



with A = 6irr]r. On time scales r » t; = m p /X 
(tj ~ 0.05 ps in our experimental setup), inertia can be 
neglected. With a linear optical force F op = —kx, this 
yields 



x(t) = ~x(t) + V2D&)r(t), 



A 



(30) 



with (r(t)r(f)) = S(t - t') and the diffusion coefficient 



k B T 



(31) 



which is known as the Einstein- Stokes equation. For 
a constant temperature equal to the room tem- 
perature of T = (294 ± 2) K, a particle diame- 
ter of 2r = (1.002 ± 0.043) pm, and a viscosity of 
r] =(1.00±0.02) -lO" 3 Nsm" 2 , the Einstein- Stokes equa- 
tion predicts a diffusion constant of 



= (0.43 ± 0.03) (pm)V 



(32) 



A more realistic treatment goes beyond Stokes' law and 
also considers the momentum that is transferred from 
the particle to the fluid. Ref. [25[ gives a derivation of 
Eqs. (|25l) . (|2"6")l . and (|27p for a macroscopic sphere in an 
incompressible fluid from linearized stochastic hydrody- 
namic equations. Thereby the Fourier transform of the 
memory kernel is computed as 



1 + 



lUJpfl 



97/ 



(33) 



Here, pf denotes the mass density of the fluid. This 
corresponds to a friction force 



F b (t) 



6Trrjrx(t) —x(t) 



6r 2 \/npfri 



dt 



(34) 



where mf — (An/3)pfr 3 is the mass of the displaced fluid. 
The correlation of the thermal fluctuations shows a neg- 
ative algebraically decaying tail 

(F th (t)F th (t')) - -ir 2 k B T^¥pJ^\t - i'|- 3/2 (35) 

for \t — i'|~ 3 / 2 > 0. The memory term in the friction 
force (f34j) and the long correlations of the thermal fluc- 
tuations clearly obliterate the Markov property of the 
process. However, it is possible that a finite Markov- 
Einstein (ME) time scale tme exists, which means that 
the process can approximately be described by a Markov 
process for times larger than tme- If such a time scale is 
found, one can try to find a model in form of Eq. (|30[) 
with an effective drag coefficient A and diffusion coeffi- 
cient Z)( 2 ) with our data analysis method. 

To obtain a rough estimation of the ME time scale, 
we compare the mean squared displacement MSD(t) = 
((Ax(t)) 2 ) of the full hydrodynamic model, which was 
computed by Clercx [2f|, with the MSD of the Markov 
model, Eq. Fig. [S] shows a plot of the MSDs 

for both models with the parameters according to our 
experiment. In the overdamped Markov model (red 
dashed curve) , there is only one characteristic time scale 
Tfe = X/k ~ 10~ 4 s depending on the stiffness k of the op- 
tical trap. Below this time scale there is the so called dif- 
fusive regime in which the MSD grows linear with time, 
MSD M (t < T fe ) = 2D^ 2 h. Above this time scale the 
MSD saturates to the constant value 2k B T/k. 




FIG. 6. (Color online) Mean squared displacement against 
time according to the full hydrodynamic model (blue solid 
line) and the overdamped Markov model (red dashed line) 
as a double logarithmic plot. The inset shows the relative 
deviation, Eq. (|36|) . between both curves. 



In the full hydrodynamic model the smallest charac- 
teristic time scale is the inertia time scale tj ~ 10 -8 s. 



8 



For times smaller than 77 the MSD grows quadratically 
with time, MSD H (t < 77) = (k B T/m*)t 2 , where m* = 
TO p + TO//2. This is the so called ballistic regime [13] ■ 
Above the ballistic regime it takes about four decades 
until the MSD approaches the diffusive regime and co- 
incides with the overdampcd Markov model. These four 
decades are influenced by hydrodynamic memory effects. 
The inset of Fig. [5] shows the relative deviation 



d{t) 



\MSD H (t)-MSD M (t)\ 
MSD H (t) 



(36) 



between the MSDs of the two models. According to this, 
the influence of the hydrodynamic memory effects is still 
present at times t ~ Tfc. For times larger than 0.4 ms, 
the relative deviation becomes smaller than one percent. 
Therefore, we expect a ME time scale of this order of 
magnitude. A direct test of our data for Markovianity 
is presented in Sec. IIVDI yielding a ME time scale of 
tme ~ 0.5 ms which is in good agreement to the discus- 
sion above. 



C. Data preparation 

Fig. [7] (c) shows the first twelve seconds of the x com- 
ponent of the original data set. One notices that the 
mean value of the time series fluctuates over time, i.e., 
the time series is not stationary. These fluctuations of 
the mean value lead to large values in the power spectral 
density (PSD) of the time series at low wave numbers k 
as one can see in Fig. [7] (a) , which shows the PSD aver- 
aged over the 150 measured trajectories. To enhance the 
stationarity of the time series, the corresponding Fourier 
coefficients for each of the 150 trajectories are reduced 
such that they fit in the Lorentz-like form of the spec- 
trum. The zero wave number coefficients are set to zero 
to center time series around zero. This kind of high pass 
filtering leads to the averaged power spectral density de- 
picted in Fig. [3(b). The corresponding time series after 
this preprocessing is shown in Fig. [3(d). The same type 
of preprocessing is applied to the y component of the 
trajectories. 



D. Preinvestigations 

As a first preinvestigation, we take a look at the 
Markov property. A necessary condition for a process 
to be Markovian on a specific time scale r is the validity 
of the Chapman-Kolmogorov equation (CKE) [l[: 



P2r{x'\x)= / dx"p T (x'\x")p T {x"\x) . 



(37) 



where p T (x'\x) := p(x',t + r\x,t) denotes the transition 
probability density function of the process. We further 
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FIG. 7. Preprocessing of the data set. Upper panels: Av- 
eraged power spectral density (APSD) of the x coordinate 
of the particle motion before (a) and after (b) preprocessing 
in arbitrary units. The insets show enlargements of the low 
wave number part, where the filtering occurs. Lower panels: 
Excerpt of the corresponding time series before (c) and after 
(d) preprocessing. 



define 

p™{x'\x) := j dx"p T/2 (x'\x")p T/2 (x"\x) , (38a) 

f;-(x',x):=p T -(x'\x)f(x), (38b) 

To test our data for Markovianity on a time scale r, we 
compare estimates for the joint PDFs /2 T and /^. K (cf. 
Ref. [6]). Fig. M shows the corresponding contour plots 
for t — t s and r = 2t s for the two components of the 
process. For t — t s (upper panels), one can see clear de- 
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viations between the contour lines which indicate that the 
process is not Markovian on this time scale. For r = 2r s 
the deviations vanish. This leads to the conclusion that 
the process has a ME time scale tme ~ 2r s . Therefore we 
only include conditional moments with time increments 
t > 2t s into the minimization of Eq. (|19p . 

As a next step, we take a look at the autocorrelation 
functions (ACF) 



C x (t) 



Cy(r) 



(X(t)X(t + r)) t 

{{X{t)?)t 
(Y(t)Y(t + r)) t 



(39a) 
(39b) 



to get an impression of the typical time scales of the 
system and to decide whether the sampling interval r s is 
sufficiently small for a reliable KM analysis. The ACFs 
are shown in Fig. [5] To evaluate the relaxation time as 
a typical time scale of the system, we fit an exponential, 
e _CT , to the first points of ACF and take t# « c _1 as 
a rough estimate. According to this estimate, tr ~ 2r s . 
This means that the relaxation time of the process is 
approximately equal to the ME time scale. Therefore, 
according to the discussion in Sec. IIII C[ the KM analysis 
should be possible. 

As a last step, we check whether it is possible or not 
to regard the x and y components of the particle mo- 
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FIG. 9. Autocorrelation functions for the x (top) and y (bot- 
tom) component of the particle motion. To obtain a rough 
estimation of the relaxation time tr, we fit an exponential 
exp(— cr) to the first points of the correlation function. This 
leads to the estimate tr ~ c _1 . 
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FIG. 8. (Color online) Test of the validity of the CKE for a 
time scale equal to the sampling interval (upper panels) and 
twice the sampling interval (lower panels) of the x compo- 
nent (left panels) and y component (right panels) of the data 
set. The plots show contour lines of estimated joint PDFs, 
according to Eqs. (|38[) . In the lower panels the contour lines 
match well in contrast to the upper panels indicating a ME 
time scale of tme ~ 2r s . 
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FIG. 10. Contour plot of the x component (left) and the y 
component (right) of the measured drift vector field. As one 
can see, the x component does not significantly depend on y 
and vice versa. 



tion as two independent processes. To this end, we 
measure the finite time drift vector field D^(x 7 y,T) — 
(Di 1 \x,y : T),D ( y 1) (x,y,T)) T with 

DP(x, y, t) = -(X(t + r) - X(t)\X(t) = x; Y(t) = y) 

T 

D<p{x, y, t) = -(Y(t + t) - Y(t)\X(t) - x- Y(t) = y) 

for t = t s . We assume that the qualitative form of the 
drift vector field is not affected that much from finite 
time effects and deviations from the Markov property. 
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In Fig. [10] we show contour plots of (x, y, r) (left) 
and Df' (x, y, r) (right). In a region with a radius of 
about 0.04 um around the origin of the coordinate sys- 
tem, where most of the data points are, the contour lines 
of (x, y, t) and Dy 1 ^ (x, y, r) are approximately linear 
and parallel to the y axis and x axis, respectively. Hence, 

D^(x,y,r)=Di 1 \x,r), (40a) 
D^(x,y,r)=D^(y,T). (40b) 

Therefore, we treat the two components of the particle 
motion as independent processes. 

E. Analysis results 

To estimate the drift and diffusion coefficients for each 
of the two processes, we first make a parametric ansatz 
in form of an OU process, Eqs. (|13p . with optimization 
parameters 7 and a. Our experience from synthetic time 
series data is that the best results are achieved if one 
includes the finite time coefficients of the n smallest time 
increments Tj = ir s , such that t„ is between one and two 
relaxation times. Since we have a finite ME time scale of 
2t s in our case, we include the conditional moments with 
t = ir s , i = 2,3,4 into the least squares potential, Eq. 

CE3). 



TABLE II. Results for the optimization parameters 7 and a 
for the x and y components of the process together with the 
error estimates umcep obtained by the MCEP method. 





7 [s- 1 ] 


a [(umfs- 1 ] 


x comp. 


result 


2004.6 


0.32910 


(J MCEP 


4.2 


0.00064 


y comp. 


result 


2212.4 


0.32240 


OMCEP 


4.8 


0.00067 



Tabic HTJ shows the results of the optimization as well as 
the estimated errors by the MCEP method for the x and 
y components of the process. A graphical representation 
is depicted in Figures ITT1 and [T2l respectively. A possible 
spatial dependence of the temperature due to the heat- 
ing of the laser cannot be resolved on the basis of the 
experimental data. If one includes a quadratic term in 
the diffusion ansatz as in Eq. (I24b[) . the error estimate 
for the parameter /3 is of the same size as the estimated 
value. However, our model with linear drift and constant 
diffusion describes the process very well, as we will see in 
Sec. ITVGl 

The diffusion coefficients deviate by a factor of approx- 
imately 1.3 from the result that was expected from the 
Stokes-Einstein equation (s. Sec. IIVB[) . To understand 
this deviation, we also measure the diffusion coefficients 
of different freely diffusing particles. The results are pre- 
sented in the following section. The stiffness of the trap 
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FIG. 11. Result of the optimization for drift (top) and diffu- 
sion (bottom) coefficients of the x coordinate of the particle 
motion. The symbols with the error bars are the estimated 
finite time coefficients for the smallest available time incre- 
ment above the ME time scale, i.e. r = 2r s . The optimized 
coefficients are represented by the solid lines. 

can nevertheless be calculated by 

k = k B T-. (41) 
a 

If we assume a temperature of (294±2) K, we obtain 

pN 

k x = (24.72 ±0.27)— , (42) 
um 

pN 

k y = (27.85 ±0.31)— . (43) 



F. Diffusion coefficient of freely diffusing particles 

In order to understand the deviations from the 
Einstein- Stokes equation, we also determine the diffu- 
sion coefficients for freely diffusing particles. In the same 
manner as described in Sec. IIV A[ we measure the posi- 
tions of seven different particles that are not trapped by 
optical tweezers. For each particle, one time series with 
approximately 10 4 time steps at a sampling frequency of 
3873 Hz is measured. From this data we determine the 
mean squared displacement (MSD) for which the relation 

MSD(t) := ((x(t) - x(0)) 2 ) = 2D (2 h (44) 

holds. Fig. rj~3"lshows the obtained MSDs for the x compo- 
nents of the seven particles. The diffusion coefficients can 



11 




-0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 
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FIG. 12. Result of the optimization for drift (top) and diffu- 
sion (bottom) coefficients of the y coordinate of the particle 
motion. The representation is analogous to Fig. [11] 

be determined by linear fits. If one averages the deter- 
mined coefficients over all particles and both coordinate 
directions, one obtains 

D {2) = (0.44 ± 0.06) (Lim) V 1 , (45) 

which is in good agreement to the expected result ac- 
cording to the Einstein-Stokes equation, Eq. The 
reason for the higher standard deviation of 0.06 (u.m) 2 s _1 
might be that the fluctuations among the particle radii 
are larger than indicated by the manufacturer. However, 
the low diffusion coefficient found in the optical trapping 
experiment is inside the range of fluctuations of diffusion 
coefficients among different particles. 

G. Comparison between model and data 

In this section we compare our estimated models to 
the experimental data. At first we compare the single 
point PDFs. The PDFs of the experimental trajectories 
are estimated via a standard kernel method using the 
Epanechnikov kernel ([23]) and a bandwidth according to 
Eq. ([15)1 . The PDFs according to the OU model are 
given by 




Fig. [H] shows the estimated PDFs from the experi- 
mental data in comparison to the ones from the model 




0.0005 0.001 0.0015 0.002 0.0025 

t[s] 



FIG. 13. Mean squared displacement of the x coordinate as 
a function of time for seven freely diffusing particles. The 
diffusion coefficient can be determined by a linear fit for each 
particle. The solid straight line shows the MSD that is ex- 
pected according to the Einstein-Stokes equation, Eq. (|32p . 

time series for the two processes. In both cases the PDFs 
agree very well. 
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FIG. 14. Comparison between the PDFs of the experimental 
data set (kernel density estimate) and our model (Eq. (|46I0 . 
for the x component (points and solid line, left vertical axis) 
and y component (squares and dashed line, right vertical axis) 
of the particle motion. 

The second quantity we compare is the ACF, Eqs. 
([39]) . Fig. [15] shows the ACF of the experimental time 
series for the x and y components of the particle motion, 
together with ACFs of our model that are given by 

C(r) = cxp(- 7 r) . (47) 

As one can also see in Fig. [S] the experimental ACFs 
are not exactly exponentially shaped. The reason for 
this is discussed in Sec. [V] However, one can see that 
the relaxation times of the models and the experimental 
data sets coincide very well. 

As a next step, we compare the first and second finite 
time conditional moments. They are depicted in Figs. 
1161 and [T71 respectively for the x component. The corre- 
sponding figures for the y component are not shown, but 
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FIG. 15. Comparison between the measured autocorrelation 
functions of the data set and the ACF of the OU model, Eq. 
(|47|l for the x component (points and solid line, left vertical 
axis) and y component (squares and dashed line, right vertical 
axis) of the particle motion. 



FIG. 17. (Color online) Comparison between the second finite 
time conditional moment of the experimental data set (red 
solid lines) and the model reconstruction (blue dashed lines) 
for the x component of the particle motion. 
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FIG. 16. (Color online) Comparison between the first finite 
time conditional moment of the experimental data set (red 
solid lines) and the model reconstruction (blue dashed lines) 
for the x component of the particle motion. 



FIG. 18. (Color online) Comparison between the fourth finite 
time KM coefficient of the data set (red solid lines) and the 
model reconstruction (blue dashed lines) for the x component 
of the particle motion. 



are qualitatively equal, 
model are given by 



The conditional moments of our 



MW(i) = -x(l-e" 7T ) 
M^(x) =x 2 (l-e-^) 2 



-(1 

7 



-2 7 t\ 



(48a) 
(48b) 



Apparently, our model fits very well to the data. Signif- 
icant deviations do only occur for the smallest r = r s , 
which is below the ME time scale and was not included 
in the optimization. Here the slope (with respect to x) 
of the first moment and the second moment (for all x) of 
the data are smaller than in the model. An inclusion of 
time increments below the ME time scale would therefore 
lead to underestimated drift and diffusion coefficients. 



Instead of the fourth conditional moments, we compare 
the fourth finite time KM coefficients which reads for our 
model 



D wr x ) =_L 

T y ' 24r 



(1 



\4 4 
Z T ) X 



6— 



1 - 2z r 



2zl- 



■3(2 



Here we have used the abbreviation 



-•JT 



(49) 



Fig. 



IT51 shows the estimated fourth KM coefficient together 
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with Eq. (|49|) . As in the two previous plots, significant 
deviations are only visible for the smallest t — t s , where 
the coefficient estimated from data is smaller than in the 
model. However, one can see that D^\x) vanishes as r 
approaches zero. This is the requirement for the Pawula 
theorem that guarantees that also the third and all higher 
KM coefficients vanish [l|. 

V. CONCLUSION 

In the present article we have demonstrated that fi- 
nite time effects can lead to significant quantitative and 
also qualitative errors in the KM analysis. A previously 
published method [l2j allows in principle to correct for 
errors caused by finite time effects. But even the ap- 
plication of this method bears the danger of misinter- 
pretation of the achieved results if the sampling interval 
of an experimental time series data set approaches the 
limit of statistical independence discussed in [10| . In or- 
der to avoid these misinterpretations, we have extended 
this method by the calculation of error estimates for the 
determined model parameters. These error estimates al- 
low for a more honest assessment of the validity of the 
obtained results. Since the estimated errors only take 
into account errors caused by finite time effects and the 
finite amount of available data, they should be regarded 
as lower bounds for the true model errors, which can also 
be influenced by other effects. 

As long as analytic solutions of the AFPE for the se- 
lected parametrization of the KM coefficients are known, 
the computational effort is very low. The applications 
presented in this articles have a computation time of less 
than five seconds on a usual desktop computer. If the 
AFPE has to be solved numerically, the computational 
effort increases dramatically. The MCEP method has not 
been tested for those cases. 

To apply our method to real-world stochastic data, we 
have also conducted an experiment where trajectories of 



Brownian particles trapped by an optical tweezers sys- 
tem were measured. We find a ME time scale which is 
approximately equal to the relaxation time of the pro- 
cess. This can be explained by hydrodynamic memory 
effects that are still present on this time scale. The large 
ME time has the consequence that even if the trajectories 
of the particle were measured with a higher sampling fre- 
quency, finite time effects could not be reduced, because 
time increments below the ME time scale must not be 
regarded in the KM analysis. 

On time scales above the ME time scale, the process 
can almost perfectly be reconstructed by an OU process 
according to the classical overdamped Markov model of 
Brownian motion. The data quality does not allow to de- 
tect deviations from linearity in the drift term or a spatial 
dependence of the diffusion one could expect because the 
laser heats up the fluid. 

The size of the measured diffusion coefficient is about 
1.3 times smaller than the diffusion predicted by the 
Einstein-Stokes equation. For comparison, we also mea- 
sured trajectories of different freely diffusing particles 
with our experimental setup. Averaged over all parti- 
cles, the Einstein- Stokes equation was found to be valid. 
The fluctuations of diffusion constants among different 
particles, which can probably be traced back to fluctu- 
ations among particle radii, are large enough to explain 
the low diffusion found for the trapped particle. 
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